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Two approaches are used to extend ENO schemes to treat conservation laws with stiff source 
terms. One approach is the application of Strang’s time-splitting method. Here the basic ENO 
scheme and Harten’s modification using subcell resolution, ENO/SR scheme, are extended this 
way. The other approach is a direct method and a modification of ENO/SR. Here the technique 
of ENO reconstruction with subcell resolution is used to locate the discontinuity within a cell and 
the time evolution is then accomplished by solving the differential equation along characteristics 
locally and advancing in the characteristic direction. This scheme is denoted ENO/SRCD. All 
the schemes are tested on the equation of LeVeque and Yee (NASA TM 100075, 1988) model- 
ing reacting flow problems. Numerical results show that these schemes handle this intriguing 
model problem very well, especially with ENO/SRCD which produces perfect resolution at the 
discontinuity. 



1. INTRODUCTION 

Recently, in studying numerical methods for reacting flow problems, LeVeque and Yee (ref.5) 
considered certain fundamental questions concerning the extension of current finite difference 
techniques developed for non-reacting flows to reacting flows. Namely, can one: (i) develop stable 
methods, (ii) obtain “high resolution” results with sharp discontinuities and second order accuracy 
in smooth regions, and (iii) obtain the correct jumps at the correct locations? They introduced 
and studied the following one-dimensional scalar conservation law with parameter-dependent 
source term 


u t + u x = ip(u), (1) 

H u) = -//u(u-~)(u-l ), (2) 

where /z is a parameter. This equation becomes stiff when the parameter /z is large. Although 
this linear advection equation with a source term represents only a simple model of reacting flow 
problems, by studying its numerical solutions one would encounter some of the difficulties sure 
to occur in solving more realistic models. 
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In their study, two different approaches were used to construct second order accurate nu- 
merical methods. One approach was to use a modification of MacCormack’s predictor-corrector 
method for conservation laws (ref.6), together with two TVD-like versions with appropriate lim- 
iters (ref.8,9). The other approach was based on the second order accurate Strang’s time-splitting 

method (ref.7). Their numerical tests revealed that stable and second order schemes can be de- 
vised by using either of these approaches. However, in studying the ability of these methods in 
dealing with propagating discontinuities, it was reported that for a certain reasonably fixed mesh 
and for the very stiff case, all the methods produced solutions that look reasonable and yet are 
completely wrong, because the discontinuities are in the wrong locations. Their investigation 
pointed out that the main difficulty is the smearing of the discontinuity in the spatial direction, 
which in turn introduced a nonequilibrium state into the calculation. To avoid this difficulty, 
it will be necessary to increase the resolution near the discontinuity, at least for the purpose of 
evaluating tp[u). 

The purpose of this paper is to show that numerical methods can be devised to overcome 
the above mentioned difficulties. We will demonstrate this by extending ENO schemes to treat 
conservation laws with source terms. We will construct numerical schemes which, when applied 
to Eq.(l), will produce stable solutions with excellent resolutions at the correct locations of 
discontinuities. We choose to describe the extensions for the following equation 

Uf\- au x = xi){u), a > 0, (3) 

where the source term ip(u) arises from the chemistry of the reacting species. It can be handled 
similarly for a < 0. 

The basic ENO scheme (ref. 4) depends on an essentially non-oscillatory reconstruction pro- 
cedure. Harten (ref.3) recently modified this procedure using subcell resolution. The notion of 
subcell resolution is based on the observation that unlike point values, cell averages of a discontin- 
uous piecewise-smooth function contain information about the exact location of the discontinuity 
within the cell. Using this observation in his study of conservation laws, Harten (ref.3) obtained 
a modification of the basic ENO scheme, which he denoted ENO/SR, achieving significant im- 
provement in the resolution of contact discontinuities. Basically, when good approximations to 
the locations of discontinuities inside the cells are obtained, it is then possible to have good 
reconstruction of the solution at each time step. Here we will also demonstrate that when the 
information on the location of the discontinuity is used in treating the source term, significant 
improvement in numerical results can be obtained. 

One approach that we use to extend ENO and ENO/SR is the application of Strang’s time- 
splitting method (ref. 7), in which one alternates between solving the conservation law without the 
source term and the ordinary differential equation modeling the chemistry. The other approach is 
a direct method and a modification of ENO/SR. Here we use the technique of ENO reconstruction 
with subcell resolution to locate the discontinuity within a cell and then accomplish the time 
evolution by solving the differential equation along characteristics locally and advancing in the 
characteristic direction. Since the subcell resolution and the characteristic direction are essential 
in the design of this scheme, it is denoted ENO/SRCD. In ref.l, ENO/SRCD was originally 
implemented using the time-splitting method. 

In section 2, we will first describe the construction of ENO/SRCD and then the extensions 
of ENO and ENO/SR schemes for Eq.(3). In section 3, we present the numerical results obtained 
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from using these schemes on the model problem of LeVeque and Yee (ref.5). A conclusion will be 
given in section 4. 


2. CONSTRUCTION OF THE SCHEMES 

We observe that along the characteristic x = xo + at, the solution to (3) evolves according 
to the ODE 

— u(x 0 + at,t) = xl)(u(x 0 + at,t)), (4) 

at 

with initial data u(x 0 ,0). In the scheme ENO/SRCD, this equation will be solved approximately 
from the time step t n to t n+ 1 . At t n , suppose that we have obtained the numerical solution 
v n = {vy }, where vj represents an approximation to u” , the cell average of u over [zy_i , Zy + i] 
at t n . Then, to obtain v n+1 , we use the following steps: 

1. Obtain a reconstruction, R(x;v n ), of the solution from the given values v n . 

2. Locate the discontinuity, if any, within each cell using the subcell resolution technique and 
modify the reconstruction R(x',v n ) to obtain R(x-,v n ). 

3. Advance R(x;v n ) via the ODE (4) along the characteristics to the t n +i level and then take 
cell averages to complete v n+1 . 

In ENO/SRCD, steps 1 and 2 will follow the ENO reconstruction procedure with subcell 
resolution of Harten (ref.3). The reconstructed solution function R(x;v n ) here is a piecewise 
quadratic polynomial obtained by using the primitive function approach. For the sake of com- 
pleteness, we will describe in straightforward terms the procedures used. For more details and 
general discussions on reconstruction with subcell resolution, see ref.3. 

Step 1. ENO Reconstruction 

Over each cell [zy_i , Zy + i ], choose i such that j — 2 < i < j and minimizes the following: 
K + 2 - 2u? + i + Vi I = min{ |v£ +2 - 2 v£ +1 + wj*| : k=j- 2 ,j - 1 ,j}. 

Let Rj(x‘,v n ) denote the reconstructed quadratic polynomial over this cell. Then 

Rj(x- v n ) = ay + Sj (x - Xj) + ^ cy (x - zy) 2 , (5) 

where 

c i = i v ?+2 -2< + i + vf)/(Az) 2 , 

s j = K+i ~ v?)/Ax + (j - i - i) cy Ax, 

a j = fy — Cy (Ax) 2 /24. 

(6) 


Step 2. Modification by Subcell Resolution 
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To detect a discontinuity in a cell [xy_ i,xy + i], we define 



In our numerical experiment, the following criterion is used. If 

M > \ s 3-il M > l 5 i+il» and Fj(Xj_l) Fj(x j+ l) < 0, (8) 

we consider that there is a discontinuity at 0y in this cell satisfying 

= (9) 


The location dj can then be approximated by using any standard root-finding method. We simply 
use the bisection method in our experiment. 


Now, if there is a discontinuity inside the cell at dj, a 
used , where 



Rj-i{x\v n ), 

R j+1 (x-,v n ), 



modified reconstruction Rj[x‘,v n ) 

< x < 0j, 

< X < X,- , i . 

j ' 2 


is 


Otherwise, we use 

Rj(x]v n ) = Rj{x\v n ). 


Step 3. Time Evolution and Cell Averaging 

Consider the case aAt < Ax and that there exists a discontinuity at 0 j inside the cell 

\Xj_L,X j+ x} With 

Oj < x j+ i - a At, 

as shown in Figure 1. At the t n level, Rj- i(x;u n ) is used as the solution to the left of 0j and 
Rj+ i(x; v n ) to the right of 0j. Since the solution to Eq.(3) evolves according to the ODE (4) along 
characteristics, we can obtain approximate solution values at the t n+i level by solving Eq.(4) and 
advancing in the characteristic direction. If w(x,t n+ i) denotes the solution of Eq.(4) at £ n+1 over 
\xj_ i , Xj + a], obtained by using the initial values Rj-i (x; v n ) on [xy_i —a At, 0j ) and Rj+i (x; u n ) 
on (0j,Xj + i — a At], the numerical solution Uy +1 is then an approximation to 


1 

Ax 



w(x,t n +i) dx. 
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^n+1 


tn 


Fig. 1. The case 0y < Xy + i — a At 

In our present implementation we use the following simple computation. Let x m and x p 
denote the midpoints in the intervals (xy_i — aAt,6j) and (0y,xy + i — a At) respectively (see 
Fig. 1). Then we compute 

w m = Rj — i (®m ! ^ ) T At (x m ; v )) , 

Wp = Rj+i (x p ; v ) + At {^p] v )) ) 

Vy +1 = [w m (9j - Xy i + a At) + w p (xy + i —0j — a At)] / Ax. 


( 10 ) 

In the above computation, tu m and w p are obtained from the Euler’s method. If the modified 
Euler method is desired, for example, one simply computes 

w* m = Rj-i(x m \ t/ n ) + At t/>(i2y_i (x m ; t> n )), 

w m = Rj-i(x m ;v n ) + — [ip(Rj-i(x m \ v n )) +ip(w* m )], 

similarly for w* and w p , and then v'J +1 by (10). Other locations of Oj and the cases with regions 
without discontinuities can be treated similarly and easily. It is quite simple to modify the above 
scheme to obtain higher order versions. 

The extension of ENO and ENO/SR schemes to treat conservation laws with source terms 
will be done by using Strang’s time-splitting method (ref.7). With respect to Eq.(3), the numerical 
solution v n+1 is computed from v n by 

-" +, =s*(y)5/(a<)^(y)»", (11) 
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where Sf(At) represents the numerical solution operator for the conservation law without the 
source term 

ut + au x — 0, a > 0 , ( 12 ) 

over a time step At, and S^(^) represents the numerical solution operator for the ordinary 
differential equation 

u t = 0(u) (13) 

over At/2. Thus the basic ENO and ENO/SR are used here as the operator Sf(At). The extended 
schemes will still be denoted ENO and ENO/SR respectively. 

The following second order version of the ENO scheme has been used in ref. 2. In our com- 
putational experiment, it produces slightly better results than other second order versions. 

ENO Scheme: 

For the operator Sf (At), we use 

( 14 ) 


where 


with 


/j'+l — “ 2 [f R (vji x j+± ^n),yj+i{xj + k ,f n )) 

d" / {vj{Xj+x , t n-|- 1 ) , Vj+1 i X j+i> ^n+l))] ) 


Ax 


+ Sj, 


Ax 


Vj+l{x j+ L,t n ) = V ? +1 - — Sj + 1, 


j+i 2 
Ax 
2~ 


Vj( x j+%>tn+i) — ^ s j AtaSj, 

Ax 
~2 


Vj+ifaj+h’tn+i) ~ v j+i o A tasj + i, 


(15) 


where the Sj's used in the above computation come from (6) in step 1 and f R (vR,VR ) denotes 
the flux at the origin in a Riemann problem with vl to the left and vr to the right. 

Now, let us describe the operator S^(At). Here we will follow steps 1 and 2 to find the 
discontinuity 9j , if any. Let us use the same notations introduced before and refer to Fig.l. Also, 
let z m and z p denote the midpoints in the intervals (xy_i,0y) and (0y,Xy + i) respectively. Then, 
for approximating cell averages and for the case Oj < Xj+t — a At, we use 


At 


S^(At) vf = v? d- 2 ^ ['l>{Rj-i{zm\v n )) [Oj - Xj_L) + ip{w m ) (6j - s y _i + a At) 
+ ^(Rj+x{z p -,v n )) i x j+$ ~ 6 i) + (w P ) {x j+ L - Oj - a At)], 


where u> m and w p are the same as in (10). Again, other situations are handled similarly. 


( 16 ) 


220 



The resulting algorithm then takes the following form: 


A f /\ f 

»r 1 = s *(7) s /( i ‘) s *(y)»'- (u) 


ENO/SR Scheme: 

The operator Sf (At) is now replaced by Harten’s second order ENO scheme with subcell 
resolution (ref.3). It is given in the form of Eq.(14) with 



-.eno A 
-fj+h + 3j+i 


where /y + 1 will be the same as in (15) and the correction term (jj + i is computed as follows. If 
the discontinuity condition (8) is not satisfied, then 


ffy+i = ^ ( u ~ x ) ~ 2 ) c i ( Ax ) 2 ’ 


v = a At/ Ax; 


else 



' [(Ax — a At) (u" — aAtsj/2) — 6y_i(xy_x,xy + i — a At)] /At, 

when F A x i+i ~a At) Fy(xy_x) > 0, 
[6y +1 (xy + i - a At,Xy + i) - a At (vy- + (Ax - a At) sy/2)]/At, 

. otherwise, 


and the expression 6y(yi,y 2 ) is used to mean 

bj{yi,y2) = / A( x ; t;n ) 

^ yi 

In the above computation, all the cy’s, sy’s, and ay’s come from (6) in step 1. The operator 
<5^ (At) will be the same as in (16) and the final algorithm also takes the same form as in (17). 


3. COMPUTATIONAL RESULTS 


We use the same fixed mesh and initial data as in the model problem of LeVeque and Yee 
(ref.5) to test the ability of the above schemes in dealing with propagating discontinuities. Thus, 
we apply all the schemes to Eq.(l) together with the initial condition 


u(x,0) = 



if x < 0.3, 

if x > 0.3. 


We take Ax = 0.02, At = 0.015, and the domain in x to be from 0 to 1. For comprison with 
ref.5, we also show the results at t = 0.3 and for the cases fj, = 1, 10, 100, and 1000. For all cases, 
ENO/SRCD produces a perfect resolution as shown in Figure 2. Figure 3 shows the computed 
results using ENO and ENO/SR schemes for /j, = 1, 10, and 100. Here one can see that the 
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results from ENO/SR are also almost perfect. For the very stiff case, /i = 1000, both ENO and 
ENO/SR fail to produce stable solutions for the above mesh in our computational experiment. 
However, when we reduce the size of At to one half of the original, i.e., At = 0.0075, and march 
40 time steps, very good result is again obtained from ENO and perfect resolution is obtained 
from ENO/SR as shown in Figures 4a and 4b respectively. We understand that reducing At 
means the reduction of stiffness of the problem. The difficulty arises from the fact that in both 

UNO 

ENO and ENO/SR schemes, the computation of the numerical flux / - + 1 still produces “large ” 
error in the spatial direction. 

The computational results obtained here compare favorably to those in LeVeque and Yee 
(ref.5). 


ENO/SRCD 



Fig. 2. Numerical results at t = 0.3 using ENO/SRCD scheme with 
discontinuous initial data and \x — 1, 10, 100, 1000. 

Ax = 0.02, At = 0.015, : true solution, • • •: computed solution 
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ENO 


ENO/SR 




Fig. 3. Numerical results at t = 0.3 using ENO (first column) and ENO/SR 
(second column) schemes with discontinuous initial data for 
fi = 1 (first row), fi = 10 (second row), and n = 100 (third row). 

Ax = 0.02, At = 0.015, — : true solution, • • •: computed solution 
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ENO 



Fig. 4a. Numerical results at t = 0.3 using ENO scheme with 
discontinuous initial data and (i — 1000. 

Ax = 0.02, At = 0.0075, : true solution, • • •: computed solution 


ENO/SR 



Fig. 4b. Numerical results at t = 0.3 using ENO/SR scheme with 
discontinuous initial data and /r = 1000. 

Ax = 0.02, At = 0.0075, — : true solution, • • •: computed solution 
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4. CONCLUSIONS 


We have extended the basic ENO and Harten’s ENO/SR schemes to treat conservation laws 
with source terms. Two approaches are used. One is to apply Strang’s time-splitting method, 
in which one alternates between solving the conservation law without the source term and the 
ordinary differential equation modeling the chemistry. The other is a modification of ENO/SR 
and a direct method, which uses the technique of ENO reconstruction with subcell resolution to 
locate the discontinuity within a cell and then accomplishes the time evolution by solving the 
differential equation along characteristics locally and advancing in the characteristic direction. 
We call this scheme ENO/SRCD. All the schemes are tested on the equation of LeVeque and Yee 
(ref.5) modeling reacting flow problems. The ENO/SRCD scheme produces perfect resolution at 
the propagating discontinuity. The extensions of basic ENO and ENO/SR via time-splitting also 
perform very well, especially with ENO/SR showing almost perfect results, except for the very 
stiff case where some adjustment in the time step-size is needed. 
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